library(readstata13)
library(dplyr)
library(grf)
library(xtable)
library(sandwich)
library(PLCE)

## Flip to "TRUE" to run again ----
runall<-FALSE


## Load data ----
Data1<-read.dta13("mattesweeksAJPS.dta")
options(device="quartz")
Data2<- Data1 %>% subset(.,6-hddv1>0)
y0<-((6-Data2$hddv1)-1)/4*100
vars.control<-c("hawk","intl","trust","voted16","birthyr","gender","educ","pid7","ideo5","newsint","pew_religimp")
X<-Data2[,vars.control]
set.seed(1);X.shuff<-apply(X,2,sample)
treat1<-Data2$hawk_t==1 & Data2$rapproche_t==1
treat2<-Data2$hawk_t==1 & Data2$rapproche_t==2
treat3<-Data2$hawk_t==2 & Data2$rapproche_t==1
treat4<-Data2$hawk_t==2 & Data2$rapproche_t==2

sub1<-(Data2$hawk_t)==1
sub2<-(Data2$hawk_t)==2

  if(runall){
    set.seed(1);h1<-plce((y0>50)[sub1],treat=treat1[sub1],X=(X)[sub1,],var.type="HC0",num.fit = 1000,printevery = 25)
    set.seed(1);h2<-plce((y0>50)[sub2],treat=treat3[sub2],X=X[sub2,],var.type="HC0",num.fit = 1000,printevery = 25)
    
    output<-list("h1"=h1,"h2"=h2)
  save(output,file="mattesHOE.Rda")
  } 


load("mattesHOE.Rda")

lm1.0<-(lm((y0>50)~treat1,sub=Data2$hawk_t==1))
lm1.1<-(lm((y0>50)~treat1+as.matrix(X),sub=Data2$hawk_t==1))


lm2.0<-(lm((y0>50)~treat3,sub=Data2$hawk_t==2))
lm2.1<-(lm((y0>50)~treat3+as.matrix(X),sub=Data2$hawk_t==2))

points1<-c(output$h1$point,lm1.0$coef[2],lm1.1$coef[2])
points2<-c(output$h2$point,lm2.0$coef[2],lm2.1$coef[2])

se.hc<-function(x) vcovHC(x,"HC0")[2,2]^.5
ses1<-c(output$h1$se,se.hc(lm1.0),se.hc(lm1.1))
ses2<-c(output$h2$se,se.hc(lm2.0),se.hc(lm2.1))

uppers1<-points1+1.96*ses1
lowers1<-points1-1.96*ses1

tab1<-rbind(points1,ses1,points2,ses2)*100
rownames(tab1)<-c("Hawks","s.e.","Doves","s.e.")
colnames(tab1)<-c("PLCE","Diff-in-Mean","OLS w Covariates ")
xtable(tab1,digits=2)

